Reanalysis of No evidence that more physically attractive women have higher estradiol or progesterone.
Data and R code see https://osf.io/qd9bv
Environment setups¶
In [1]:
import os
import numpy as np
import pandas as pd
import seaborn as sns
import pymc3 as pm
import theano.tensor as tt
%pylab inline
%config InlineBackend.figure_format = 'retina'
%load_ext watermark
plt.style.use('ggplot')
In [2]:
%watermark -v -m -p numpy,pandas,seaborn,pymc3,matplotlib
Load data, display descriptive result¶
In [3]:
tbl0 = pd.read_csv('ratings_anon.csv')
tbl0.head()
Out[3]:
In [4]:
measure_label = ['progesterone', 'estradiol', 'rating', 'whr']
Nmeasure = len(measure_label)
(tbl0.groupby(['face_id'])[measure_label]
.agg(['mean', 'std'])).round(3)
Out[4]:
In [5]:
dfmean = (tbl0.groupby(['face_id'])[measure_label]
.agg(np.nanmean))
sns.pairplot(dfmean[np.isfinite(dfmean.values.mean(axis=1))], diag_kind='kde', plot_kws=dict(alpha=.25));
In [6]:
dfmean1 = (tbl0.groupby(['face_id', 'session'])[measure_label]
.agg(np.nanmean))
sns.pairplot(dfmean1[np.isfinite(dfmean1.values.mean(axis=1))], diag_kind='kde', plot_kws=dict(alpha=.05));
In [7]:
_, ax = plt.subplots(1, 1, figsize=(5, 5))
dataval = dfmean.values
for i in range(dataval.shape[1]):
dataval[~np.isfinite(dataval[:, i]), i] = np.nanmean(dataval[:, i])
empirical_corr = pd.DataFrame(np.corrcoef(dataval.T),
columns=measure_label,
index=measure_label)
sns.heatmap(empirical_corr,
cbar=False, square = True, annot=True,
linewidths=.1, cmap='viridis', ax=ax)
ax.set_title('Empirical Correlation Matrix');
Base line model on the subject average¶
In [8]:
with pm.Model() as model_mean:
# on the assumption that attractiveness, progesterone and estradiol
# are correlated, this correlation could be capture with a cov matrix
sd_dist = pm.HalfCauchy.dist(5)
packed_chol = pm.LKJCholeskyCov('chol_cov', n=Nmeasure, eta=1, sd_dist=sd_dist)
# compute the covariance matrix
chol = pm.expand_packed_triangular(Nmeasure, packed_chol, lower=True)
cov = tt.dot(chol, chol.T)
# Extract the correlation deviations etc
sd = tt.sqrt(tt.diag(cov))
corr = pm.Deterministic('corr', tt.diag(sd**-1).dot(cov.dot(tt.diag(sd**-1))))
# likelihood
mu = pm.Normal('mu', mu=dataval.mean(axis=0), sd=10., shape=Nmeasure)
obs = pm.MvNormal('obs', mu=mu, chol=chol, observed=dataval)
# samples
trace = pm.sample(1000, tune=2000, njobs=4)
In [9]:
pm.traceplot(trace, varnames=['chol_cov']);
In [10]:
_, ax = plt.subplots(1, 2, figsize=(10, 5))
sns.heatmap(empirical_corr,
cbar=False, square = True, annot=True,
linewidths=.1, cmap='viridis', ax=ax[0])
ax[0].set_title('Empirical Correlation Matrix')
posterior_corr = pd.DataFrame(trace['corr'].mean(axis=0),
columns=measure_label,
index=measure_label)
sns.heatmap(posterior_corr,
cbar=False, square = True, annot=True,
linewidths=.1, cmap='viridis', ax=ax[1])
ax[1].set_title('Latent Correlation Matrix')
plt.tight_layout();
In [11]:
plotpost = pm.plots.artists.plot_posterior_op
idxall = np.tril_indices(4, k=-1)
_, ax = plt.subplots(2, 3, figsize=(18, 6), sharex=True)
ax1 = ax.flatten()
for i in range(len(idxall[0])):
iy, ix = idxall[0][i], idxall[1][i]
trace_values = trace['corr'][:, ix, iy]
plotpost(trace_values, ax1[i], kde_plot=False, point_estimate='mean',
round_to=3, alpha_level=0.05, ref_val=0., rope=None, color='#87ceeb')
ax1[i].axvline(empirical_corr.values[ix, iy], color='r')
ax1[i].set_title('Corr('+measure_label[ix]+', '+measure_label[iy]+')')
plt.tight_layout();
Full model on observed¶
Prepare data
In [12]:
# test model on a subset of data
randsubj = tbl0['face_id'].unique()[:20]
tbl = (tbl0.loc[tbl0['face_id'].isin(randsubj).values, :]).reset_index()
tbl_ = tbl[['face_id', 'session', 'progesterone', 'estradiol']].drop_duplicates()
sns.pairplot(tbl_, hue='face_id', vars=['progesterone', 'estradiol'],
diag_kind='kde', plot_kws=dict(alpha=.25));
In [13]:
tbl = tbl0
Data format (for 1 subject)
In [14]:
tbl[tbl['face_id']==tbl['face_id'][0]]
Out[14]:
Set up latent variable (N participant * K measurment)
In [15]:
with pm.Model() as mvmodel:
# on the assumption that attractiveness, progesterone and estradiol
# are correlated, this correlation could be capture with a cov matrix
sd_dist = pm.HalfCauchy.dist(5)
packed_chol = pm.LKJCholeskyCov('chol_cov', n=Nmeasure, eta=1, sd_dist=sd_dist)
# compute the covariance matrix
chol = pm.expand_packed_triangular(Nmeasure, packed_chol, lower=True)
cov = tt.dot(chol, chol.T)
# Extract the correlation deviations etc
sd = tt.sqrt(tt.diag(cov))
corr = pm.Deterministic('corr', tt.diag(sd**-1).dot(cov.dot(tt.diag(sd**-1))))
# latent (zero mean Multivariate Normal)
# latent = pm.MvNormal('latent', mu=np.zeros((Nmeasure,)), chol=chol,
# shape=dataval.shape, testval=dataval)
# transform an uncorrelated normal:
mu0 = pm.Normal('mu0', mu=0., sd=1.,
shape=dataval.shape)
latent = tt.dot(chol, mu0.T).T
In [16]:
latent.tag.test_value.shape
Out[16]:
WHR measurement get the shortest (1 observeration per participant, with one missing).
In [17]:
subjID = np.asarray(tbl['face_id'].astype('category').cat.codes, dtype=np.int16)
# whr
whr_raw = np.asarray(tbl['whr'], dtype=np.float64)
sortID, index = np.unique(subjID, return_index=True)
whr_ = whr_raw[index]
index2 = np.isfinite(whr_)
whr = np.asarray(whr_[index2], dtype=np.float64)
whr_ID = sortID[index2]
In [18]:
_, ax = plt.subplots(1, 1, figsize=(18, 3))
ax.plot(whr);
In [19]:
# whr
with mvmodel:
whr_mu = pm.Normal('whr_mu', mu=whr.mean(), sd=1.)
whr_obs = pm.Normal('whr_obs', mu=whr_mu+latent[whr_ID, 3], sd=.01, observed=whr)
Hormone measurment, each subject with 5 observations
In [20]:
# progesterone and estradiol
tbl2 = tbl[['face_id', 'session', 'progesterone', 'estradiol']].drop_duplicates()
progesterone = np.asarray(tbl2['progesterone'], dtype=np.float64)
estradiol = np.asarray(tbl2['estradiol'], dtype=np.float64)
subjID_hm = np.asarray(tbl2['face_id'].astype('category').cat.codes, dtype=np.int16)
sessID_hm = np.asarray(tbl2['session'].astype('category').cat.codes, dtype=np.int16)
# Nsess = 5
# sessID_hm = np.repeat(np.arange(Nsess)[None,], len(sortID), 0).flatten()
In [21]:
hormone_df = pd.DataFrame(np.vstack((sessID_hm, progesterone, estradiol)).T,
columns=['session', 'progesterone', 'estradiol'])
sns.pairplot(hormone_df, hue='session', vars=['progesterone', 'estradiol'],
diag_kind='kde', plot_kws=dict(alpha=.25));
In [22]:
progest_mat = progesterone.reshape(len(np.unique(subjID_hm)), -1).T
estradi_mat = estradiol.reshape(len(np.unique(subjID_hm)), -1).T
horm_ID = subjID_hm.reshape(len(np.unique(subjID_hm)), -1).T[0]
In [23]:
_, ax = plt.subplots(2, 1, figsize=(18, 6))
ax[0].imshow(progest_mat, cmap='viridis', aspect='auto')
ax[0].axis('off')
ax[1].imshow(estradi_mat, cmap='viridis', aspect='auto')
ax[1].axis('off')
plt.tight_layout();
In [24]:
_, ax = plt.subplots(2, 1, figsize=(18, 6))
ax[0].plot(np.std(progest_mat, axis=0), label='progesterone')
ax[1].plot(np.std(estradi_mat, axis=0), label='estradiol')
plt.tight_layout();
print('Standard Deviation of hormone measurement')
In [25]:
# parameters for informative prior
hm1_emp_std = np.std(progest_mat, axis=0).mean()
hm2_emp_std = np.std(estradi_mat, axis=0).mean()
In [26]:
# progesterone and estradiol
with mvmodel:
hm1_mu = pm.Normal('mu1', mu=progest_mat.mean(), sd=100.)
hm1_std = pm.HalfCauchy('sd1', beta=5, shape=np.unique(subjID).shape, testval=hm1_emp_std)
hm1_nu = pm.HalfNormal('nu1', 10.)
hm1_obs = pm.StudentT('progest_obs',
mu=latent[horm_ID, 0] + hm1_mu,
sd=hm1_std[horm_ID],
nu=hm1_nu,
observed=progest_mat)
hm2_mu = pm.Normal('mu2', mu=estradi_mat.mean(), sd=10.)
hm2_std = pm.HalfCauchy('sd2', beta=2, shape=np.unique(subjID).shape, testval=hm2_emp_std)
hm2_nu = pm.HalfNormal('nu2', 10.)
hm2_obs = pm.StudentT('estradi_obs',
mu=latent[horm_ID, 1] + hm2_mu,
sd=hm2_std[horm_ID],
nu=hm2_nu,
observed=estradi_mat)
In [27]:
hm1_std.tag.test_value.shape
Out[27]:
In [28]:
hm1_obs.tag.test_value.shape
Out[28]:
Rating, the tallest data column
In [29]:
# ratings
Nrater = len(tbl['rater_id'].unique())
raterID = np.asarray(tbl['rater_id'].astype('category').cat.codes, dtype=np.int16)
raterSex = np.asarray(tbl['rater_sex'].astype('category').cat.codes, dtype=np.int16)
rating = np.asarray(tbl['rating'], dtype=np.float64)
# change rating to real line
rating_r = np.log(rating - 0) - np.log(8 - rating)
In [30]:
Nrating = len(np.unique(rating))
np.unique(rating)
Out[30]:
In [31]:
with mvmodel:
ratermu = pm.Normal('ratermu', mu=0., sd=1., shape=np.unique(raterID).shape)
ratersd = pm.HalfNormal('ratersd', sd=10., shape=np.unique(raterID).shape)
latentrate = latent[subjID, 2] + ratermu[raterID]
rating_obs = pm.Normal('rating_obs',
mu=latentrate,
sd=ratersd[raterID],
observed=rating_r)
In [32]:
with mvmodel:
trace = pm.sample(2000, tune=1000, njobs=4)
In [33]:
pm.traceplot(trace, varnames=['chol_cov', 'nu1', 'nu2']);
In [34]:
CONVERGENCE_TITLE = lambda: 'BFMI = {a:.2f}\nmax(R_hat) = {b:.3f}\nmin(Eff_n) = {c:.3f}'\
.format(a=bfmi, b=max_gr, c=min_effn)
def get_diags(trace):
bfmi = pm.bfmi(trace)
max_gr = max(np.max(gr_stats) for gr_stats in pm.gelman_rubin(trace).values())
min_effn = min(np.median(ef_stats) for ef_stats in pm.effective_n(trace).values())
return bfmi, max_gr, min_effn
_, ax = plt.subplots(1, 1, figsize=(7, 5))
bfmi, max_gr, min_effn = get_diags(trace)
(pm.energyplot(trace, ax=ax)
.set_title(CONVERGENCE_TITLE()));
In [86]:
ppc = pm.sample_ppc(trace, 500, model=mvmodel)
In [87]:
ppc.keys()
Out[87]:
In [88]:
_, ax = plt.subplots(1, 1, figsize=(18, 6))
ax.plot(ppc['whr_obs'].T, 'k', alpha=.01);
ax.plot(whr, 'xr');
In [89]:
_, ax = plt.subplots(1, 1, figsize=(18, 6))
ax.plot(ppc['progest_obs'].T, 'ok', alpha=.01);
ax.plot(progest_mat.T, 'xr', alpha=.5);
ax.set_ylim(-100, np.max(progest_mat)+100);
In [90]:
_, ax = plt.subplots(1, 1, figsize=(18, 6))
ax.plot(ppc['estradi_obs'].T, 'ok', alpha=.01);
ax.plot(estradi_mat.T, 'xr', alpha=.5);
ax.set_ylim(-1, np.max(estradi_mat)+1);
In [91]:
pltwidth = 250
idx = 15000 #np.random.randint(len(rating_r)-pltwidth)
pltrange = np.arange(idx, idx+pltwidth)
_, ax = plt.subplots(1, 1, figsize=(18, 6))
ax.plot(pltrange, ppc['rating_obs'].T[pltrange], 'ok', alpha=.01);
ax.plot(pltrange, rating_r.T[pltrange], 'r');
In [41]:
_, ax = plt.subplots(1, 2, figsize=(10, 5))
sns.heatmap(empirical_corr,
cbar=False, square = True, annot=True,
linewidths=.1, cmap='viridis', ax=ax[0])
ax[0].set_title('Empirical Correlation Matrix')
posterior_corr = pd.DataFrame(trace['corr'].mean(axis=0),
columns=measure_label,
index=measure_label)
sns.heatmap(posterior_corr,
cbar=False, square = True, annot=True,
linewidths=.1, cmap='viridis', ax=ax[1])
ax[1].set_title('Latent Correlation Matrix')
plt.tight_layout();
In [42]:
plotpost = pm.plots.artists.plot_posterior_op
idxall = np.tril_indices(4, k=-1)
_, ax = plt.subplots(2, 3, figsize=(18, 6), sharex=True)
ax1 = ax.flatten()
for i in range(len(idxall[0])):
iy, ix = idxall[0][i], idxall[1][i]
trace_values = trace['corr'][:, ix, iy]
plotpost(trace_values, ax1[i], kde_plot=False, point_estimate='mean',
round_to=3, alpha_level=0.05, ref_val=0., rope=None, color='#87ceeb')
ax1[i].axvline(empirical_corr.values[ix, iy], color='r')
ax1[i].set_title('Corr('+measure_label[ix]+', '+measure_label[iy]+')')
plt.tight_layout();
Alternative models¶
Using LogNormal to model the hormone measurements
In [62]:
import matplotlib.pyplot as plt
import numpy as np
import scipy.stats as st
x = np.linspace(0.0, 2000.0, 1000)
fig, ax = plt.subplots()
ax.hist(progest_mat.flatten(), 100, normed=True)
f = lambda mu, sd : st.lognorm.pdf(x, sd, scale=np.exp(mu))
plot_pdf = lambda mu, sd : ax.plot(x, f(mu, sd), label=r'$\mu$={0}, $\sigma$={1}'.format(mu, sd))
plot_pdf(np.log(progest_mat.mean()), 1)
plt.legend(loc='upper right', frameon=False)
plt.show()
In [67]:
fig, ax = plt.subplots()
x = np.linspace(0.0, 50.0, 1000)
ax.hist(estradi_mat.flatten(), 100, normed=True)
f = lambda mu, sd : st.lognorm.pdf(x, sd, scale=np.exp(mu))
plot_pdf = lambda mu, sd : ax.plot(x, f(mu, sd), label=r'$\mu$={0}, $\sigma$={1}'.format(mu, sd))
plot_pdf(np.log(estradi_mat.mean()), .6)
plt.legend(loc='upper right', frameon=False)
plt.show()
In [68]:
with pm.Model() as mvmodel2:
# on the assumption that attractiveness, progesterone and estradiol
# are correlated, this correlation could be capture with a cov matrix
sd_dist = pm.HalfCauchy.dist(5)
packed_chol = pm.LKJCholeskyCov('chol_cov', n=Nmeasure, eta=1, sd_dist=sd_dist)
# compute the covariance matrix
chol = pm.expand_packed_triangular(Nmeasure, packed_chol, lower=True)
cov = tt.dot(chol, chol.T)
# Extract the correlation deviations etc
sd = tt.sqrt(tt.diag(cov))
corr = pm.Deterministic('corr', tt.diag(sd**-1).dot(cov.dot(tt.diag(sd**-1))))
# latent (zero mean Multivariate Normal)
# latent = pm.MvNormal('latent', mu=np.zeros((Nmeasure,)), chol=chol,
# shape=dataval.shape, testval=dataval)
# transform an uncorrelated normal:
mu0 = pm.Normal('mu0', mu=0., sd=1.,
shape=dataval.shape)
latent = tt.dot(chol, mu0.T).T
# whr
whr_mu = pm.Normal('whr_mu', mu=whr.mean(), sd=1.)
whr_obs = pm.Normal('whr_obs', mu=whr_mu+latent[whr_ID, 3], sd=.01, observed=whr)
# progesterone and estradiol
hm1_mu = pm.Normal('mu1', mu=np.log(progest_mat.mean()), sd=10.)
hm1_std = pm.HalfNormal('sd1', 1., shape=np.unique(subjID).shape, testval=hm1_emp_std)
hm1_obs = pm.Lognormal('progest_obs',
mu=tt.exp(latent[horm_ID, 0] + hm1_mu),
sd=hm1_std[horm_ID],
observed=progest_mat)
hm2_mu = pm.Normal('mu2', mu=np.log(estradi_mat.mean()), sd=1.)
hm2_std = pm.HalfNormal('sd2', .5, shape=np.unique(subjID).shape, testval=hm2_emp_std)
hm2_obs = pm.Lognormal('estradi_obs',
mu=tt.exp(latent[horm_ID, 1] + hm2_mu),
sd=hm2_std[horm_ID],
observed=estradi_mat)
# ratings
ratermu = pm.Normal('ratermu', mu=0., sd=1., shape=np.unique(raterID).shape)
ratersd = pm.HalfNormal('ratersd', sd=10., shape=np.unique(raterID).shape)
latentrate = latent[subjID, 2] + ratermu[raterID]
rating_obs = pm.Normal('rating_obs',
mu=latentrate,
sd=ratersd[raterID],
observed=rating_r)
In [69]:
with mvmodel2:
trace2 = pm.sample(2000, tune=1000, njobs=4)
In [70]:
_, ax = plt.subplots(1, 1, figsize=(7, 5))
bfmi, max_gr, min_effn = get_diags(trace2)
(pm.energyplot(trace2, ax=ax)
.set_title(CONVERGENCE_TITLE()));
In [71]:
pm.traceplot(trace2, varnames=['chol_cov', 'sd1', 'sd2']);
In [73]:
import pickle
tracedump = dict(trace=trace,
trace2=trace2)
with open('trace_temp.pickle', 'wb') as f:
pickle.dump(tracedump, f)
# with open('trace_temp.pickle', 'rb') as f:
# tracedump = pickle.load(f)
# trace1 = tracedump['trace']
# trace2 = tracedump['trace2']
In [76]:
# MODEL_NAME_MAP = {
# 0: "Student-T",
# 1: "LogNormal"
# }
# comp_df = (pm.compare([trace, trace2],
# [mvmodel, mvmodel2],
# )
# .rename(index=MODEL_NAME_MAP)
# .loc[MODEL_NAME_MAP.values()])
# comp_df
In [77]:
ppc = pm.sample_ppc(trace2, 500, model=mvmodel2)
In [78]:
ppc.keys()
Out[78]:
In [79]:
_, ax = plt.subplots(1, 1, figsize=(18, 6))
ax.plot(ppc['whr_obs'].T, 'k', alpha=.01);
ax.plot(whr, 'xr');
In [80]:
_, ax = plt.subplots(1, 1, figsize=(18, 6))
ax.plot(ppc['progest_obs'].T, 'ok', alpha=.01);
ax.plot(progest_mat.T, 'xr', alpha=.5);
ax.set_ylim(-100, np.max(progest_mat)+100);
In [81]:
_, ax = plt.subplots(1, 1, figsize=(18, 6))
ax.plot(ppc['estradi_obs'].T, 'ok', alpha=.01);
ax.plot(estradi_mat.T, 'xr', alpha=.5);
ax.set_ylim(-1, np.max(estradi_mat)+1);
In [85]:
pltwidth = 250
idx = 15000 #np.random.randint(len(rating_r)-pltwidth)
pltrange = np.arange(idx, idx+pltwidth)
_, ax = plt.subplots(1, 1, figsize=(18, 6))
ax.plot(pltrange, ppc['rating_obs'].T[pltrange], 'ok', alpha=.01);
ax.plot(pltrange, rating_r.T[pltrange], 'r');
In [83]:
_, ax = plt.subplots(1, 2, figsize=(10, 5))
sns.heatmap(empirical_corr,
cbar=False, square = True, annot=True,
linewidths=.1, cmap='viridis', ax=ax[0])
ax[0].set_title('Empirical Correlation Matrix')
posterior_corr = pd.DataFrame(trace2['corr'].mean(axis=0),
columns=measure_label,
index=measure_label)
sns.heatmap(posterior_corr,
cbar=False, square = True, annot=True,
linewidths=.1, cmap='viridis', ax=ax[1])
ax[1].set_title('Latent Correlation Matrix')
plt.tight_layout();
In [84]:
plotpost = pm.plots.artists.plot_posterior_op
idxall = np.tril_indices(4, k=-1)
_, ax = plt.subplots(2, 3, figsize=(18, 6), sharex=True)
ax1 = ax.flatten()
for i in range(len(idxall[0])):
iy, ix = idxall[0][i], idxall[1][i]
trace_values = trace2['corr'][:, ix, iy]
plotpost(trace_values, ax1[i], kde_plot=False, point_estimate='mean',
round_to=3, alpha_level=0.05, ref_val=0., rope=None, color='#87ceeb')
ax1[i].axvline(empirical_corr.values[ix, iy], color='r')
ax1[i].set_title('Corr('+measure_label[ix]+', '+measure_label[iy]+')')
plt.tight_layout();
Ordered logistic for Rating
In [ ]:
class Ordered(pm.distributions.transforms.ElemwiseTransform):
name = "ordered"
def forward(self, x):
out = tt.zeros(x.shape)
out = tt.inc_subtensor(out[0], x[0])
out = tt.inc_subtensor(out[1:], tt.log(x[1:] - x[:-1]))
return out
def forward_val(self, x, point=None):
x, = pm.distributions.distribution.draw_values([x], point=point)
return self.forward(x)
def backward(self, y):
out = tt.zeros(y.shape)
out = tt.inc_subtensor(out[0], y[0])
out = tt.inc_subtensor(out[1:], tt.exp(y[1:]))
return tt.cumsum(out)
def jacobian_det(self, y):
return tt.sum(y[1:])
with pm.Model() as mvmodel3:
# on the assumption that attractiveness, progesterone and estradiol
# are correlated, this correlation could be capture with a cov matrix
sd_dist = pm.HalfCauchy.dist(5)
packed_chol = pm.LKJCholeskyCov('chol_cov', n=Nmeasure, eta=1, sd_dist=sd_dist)
# compute the covariance matrix
chol = pm.expand_packed_triangular(Nmeasure, packed_chol, lower=True)
cov = tt.dot(chol, chol.T)
# Extract the correlation deviations etc
sd = tt.sqrt(tt.diag(cov))
corr = pm.Deterministic('corr', tt.diag(sd**-1).dot(cov.dot(tt.diag(sd**-1))))
# latent (zero mean Multivariate Normal)
# latent = pm.MvNormal('latent', mu=np.zeros((Nmeasure,)), chol=chol,
# shape=dataval.shape, testval=dataval)
# transform an uncorrelated normal:
mu0 = pm.Normal('mu0', mu=0., sd=1.,
shape=dataval.shape)
latent = tt.dot(chol, mu0.T).T
# whr
whr_mu = pm.Normal('whr_mu', mu=whr.mean(), sd=1.)
whr_obs = pm.Normal('whr_obs', mu=whr_mu+latent[whr_ID, 3], sd=.01, observed=whr)
# progesterone and estradiol
hm1_mu = pm.Normal('mu1', mu=np.log(progest_mat.mean()), sd=10.)
hm1_std = pm.HalfNormal('sd1', 1., shape=np.unique(subjID).shape, testval=hm1_emp_std)
hm1_obs = pm.Lognormal('progest_obs',
mu=tt.exp(latent[horm_ID, 0] + hm1_mu),
sd=hm1_std[horm_ID],
observed=progest_mat)
hm2_mu = pm.Normal('mu2', mu=np.log(estradi_mat.mean()), sd=1.)
hm2_std = pm.HalfNormal('sd2', .5, shape=np.unique(subjID).shape, testval=hm2_emp_std)
hm2_obs = pm.Lognormal('estradi_obs',
mu=tt.exp(latent[horm_ID, 1] + hm2_mu),
sd=hm2_std[horm_ID],
observed=estradi_mat)
# ratings
ratermu = pm.Normal('ratermu', mu=0., sd=1., shape=np.unique(raterID).shape)
ratersd = pm.HalfNormal('ratersd', sd=10., shape=np.unique(raterID).shape)
latentrate = latent[subjID, 2] + ratermu[raterID]*ratersd[raterID]
# latent rating
a = pm.Normal('a', 0., 10.,
transform=Ordered(),
shape=Nrating,
testval=np.arange(Nrating) - Nrating/2)
pa = pm.math.sigmoid(tt.shape_padleft(a) - tt.shape_padright(latentrate))
p_cum = tt.concatenate([
tt.zeros_like(tt.shape_padright(pa[:, 0])),
pa,
tt.ones_like(tt.shape_padright(pa[:, 0]))
], axis=1)
p = p_cum[:, 1:] - p_cum[:, :-1]
rating_obs = pm.Categorical('rating_obs', p, observed=rating-1)